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A new computational technique based on the symbolic description utilizing kneading invariants 
is proposed and verified for explorations of dynamical and parametric chaos in a few exemplary 
systems with the Lorenz attractor. The technique allows for uncovering the stunning complexity 
I and universality of bi-parametric structures and detect their organizing centers - codimension- 

^ two T-points and separating saddles in the kneading-based scans of the iconic Lorenz equation 

00 from hydrodynamics, a normal model from mathematics, and a laser model from nonlinear 

optics. 
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1. Introduction 

A great deal of analytical and experimental studies, including modeling simulations, have been focused 
on the identification of key signatures to serve as structural invariants that would allow dynamically alike 
nonlinear systems with chaotic dynamics from diverse origins to be united into a single class. Among these 
key structures are various homoclinic and heteroclinic bifurcations of low codimensions that lie at the 
heart of the understanding of complex behaviors because of their roles of organizing centers of dynamics 
in parameterized dynamical systems. 

Dynamical systems theory has aimed to create purely abstract approaches that are further proceeded 
by creation and the development of applicable tools designed for the search and identification of such 
basic invariants for simple, Morse-Smale, systems and ones with complex chaotic dynamics. One such a 
[computationally justified] approach for studying complex dynamics capitalizes on the concept of sensitivity 
of deterministic chaos. Sensitivity of chaotic trajectories can be quantified in terms of the divergence rate 
evaluated through the largest Lyapunov characteristic exponent. The approach has been proven to work 
exceptionally well for various systems with chaotic and simple dynamics. In several low-order dissipative 
systems, like the Rossler model, the computational technique based on the largest Lyapunov characteristic 
exponent reveals that they possess common, easily recognizable patterns involving spiral structures in bi- 
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parametric p] 


anes jAfraimovich et al, 1977; Bykov 


, 1993 


; Shilnikov, 1991; Shilnikov et al. 


, 1993 


; Barrio 


et al. , 


2011b; 


G alias, 2010 . Such patterns turn out to be ubiquitously alike in both time-discrete 


Lorenz , 



2008 1 and time-continuous systems |Gaspard et al. , 1984; Barrio et al., 2009; G alias, 2010 , and they are 



easily located when the spiral patterns have regular and chaotic spiral "arms" in the systems with the 
Shilnikov saddle- focus |Shilnikov fc Shilnikov , 2007 



Application of the Lyapunov exponents technique fails, in general, to reveal fine structures embedded 
in the bi-parametric scans of Lorenz-like systems. As such it cannot deliver the desired insights into intrinsic 
bifurcations because regions of chaotic dynamics appear to be uniformly. This basically means that the 
instability of the Lorenz attractors does not vary noticeably as control parameters of the system are varied. 
This holds true too when one attempts to find the presence of characteristic spiral structures that are known 
to exist theoretically in the Lorenz-like systems |Bykov , 1993; Glendinning & Sparrow, 1986| and therefore 



could only be identified using accurate bifurcation continuation approaches | Shilnikov, 1991 ; Shilnikov et al. 



1993 . Such spirals in a bi-parametric parameter plane of the system in question are organized around the 



so-called T[erminal] -points corresponding to codimension-two or -higher heteroclinic connections between 
two or more saddle equilibria. For ^-symmetric systems with the Lorenz attractor, the degeneracy is 
reduced to a cod-2 bifurcation of a closed heteroclinic connection involving two saddle-foci and a saddle at 
the origin. T-points have been located in various models of diverse origins including electronic oscillators 
Bykov, 1998[ jFernandez- Sanchez et al. , 2002| and nonlinear optics |Forysiak et al. , 1991; Wieczorek &; 



Krauskopij|2005l , etc. 
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Figure 1. Caricature of the bifurcation unfolding of the ordinary T-point for symmetric Lorenz-like systems with a closed 
heteroclinic connection evolving both saddle-foci and the saddle. Point out that with each revolution approaching the T-point 
along the curve (1-hom O), the number of turns of the one-dimensional separatrix of the saddle, O, around the saddle-focus 
increases by one in the homoclinic loop and becomes infinite at the T-point. 



Figure [T] sketches an idea of the structure of the bifurcation unfolding near an ordinary T-point in 



the parameter plane in a symmetric system |Bykov 


, 1993 


; Glendinning &; Sparrow 




1986 ; various T- 


points configurations for other heteroclinic connections were examined in detail in |Bykov , 


1993 . Here, the 



heteroclinic connection is formed by a pair of symmetric saddle-foci and a saddle whose one-dimensional 
stable (incoming) and unstable (outgoing), respectively, separatrices merge only at the codimension-two 
T-point in the bi-parametric space of the Lorenz model. It follows from the theoretical analysis that the 
unfolding of this bifurcation includes three main curves ending at the T-point: a spiraling bifurcation curve, 
(1-hom O) corresponding to two simultaneous (due to the symmetry) homoclinic loops of the saddle, a 
curve, (1-hom Oip), corresponding to two simultaneous homoclinic loops of the saddle-foci, and another 
codimension-one bifurcation curve, (1-het Oip), corresponding to a heteroclinic correction between both 
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saddle- foci. In addition, the unfolding includes infinitely many subsidiary T-points in between 1-hom 
and 1-het 0\^- 

Despite a rather overwhelming number of studies reporting the occurrence of various spiral structures, 
there is yet unproportionately little known about construction details and generality of underlying bifur- 
cation scenarios giving rise to such patterns. In this paper we introduce a novel computational toolkit 
capitalizing on the idea of the symbolic representation for the dynamics of Lorenz-like systems that em- 
ploys the kneading invariants |Milnor fc Thurston, 1988] . We will then show how the toolkit can be used for 



detecting various structures in bi-parametric scans of such systems. It is our intention to enhance further 
the technique thus allowing for systematic studies of the stunning complexity and universality of spiral 
structures in models of diverse dynamics and origins. 

The paper is organized as follows: in Section [2] we review the homoclinic bifurcation theory and discuss 
the routes to chaos and the formation stages of the strange attractor in the Lorenz model; in Section [3] 
we introduce basics of kneading theory; in Section [4] we demonstrate of the computational technique using 
kneading invariants to reveal hidden structures of biparametric chaos in the iconic Lorenz model; in Section 
[5]we apply the technique to uncover colorfulness of the bifurcation diagrams of the Shimizu-Morioka model. 
Finally, In Section [6] the proposed technique will be tested on a 6D model of the optically pumped, far 



infrared red three-level laser |Moloney et al , 1989; Forysiak et al, 1991] to confirm the universality of the 



patters produced by the deterministic chaos in the Lorenz like systems. 
2. Homoclinic bifurcations in systems with the Lorenz attractor 

The strange chaotic attractor in the Lorenz equation from hydrodynamics has become a de-facto proof of 
deterministic chaos. The butterfly-shaped image of the iconic Lorenz attractor, shown in Fig.[5j has become 
the trademark of Chaos theory and Dynamical Systems. This theory elaborates on complex trajectory 
behaviors in nonlinear systems from mathematics, physics, life sciences, finance, etc. Universality of the 
methods along with bifurcation tools has made them spread wide and deep across all other disciplines of 
the modern science. 



The Lorenz equation Lorenz, 1963 1 is a system of three differential equations: 



x — —cr(x — y), y — r x — y — xz, z — —bz + xy, (1) 

with three positive bifurcation parameters: a being the Prandtl number quantifying the viscosity of the 
fluid, b being a positive constant of magnitude order 1 which originates from the nonlinearity of the 
Boussinesq equation, and r being a Reynolds number that characterizes the fluid dynamics. Notice that 
Eqs. ([!]) are ^-symmetric, i.e. <H> (— x, — y^z) |Sparrow, 1982 . 



2.1. Uni-parametric cut through the Lorenz equation 



; \ unstable limit cycle 




Figure 2. Sketch of the uni-parametric bifurcation diagram for the Lorenz equation at a — 10 and b = 8/3: plotted are the 
coordinates, \x\, of the limit trajectories (equilibria, periodic and homoclinic orbits) against the bifurcation parameter r. 



An exploration of primary bifurcations in the Lorenz equation begins with a single-parameter exam- 
ination of the dynamics, where r serves as the bifurcation parameter increasing from laminar to weekly 
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turbulent magnitudes around 25, while the two other parameters are set to the original Saltzman values: 
a = 10 (for the water, and a = 1 for the air) and 6 = 8/3. This would give a uni-parametric cut through 



the Lorenz equation that was originally explored in these independent studies |Afraimovich et al, 1977 
Kaplan fc Yorke} (1979) (see Fig. §): 



For r < 1 the only equilibrium state O(0, 0, 0) at the origin is a global attractor in the 3D phase space of 
the Lorenz equation. 

This equilibrium state undergoes a pitch- fork bifurcation at r = rp = 1, and for r > 1 becomes a saddle 
so that the stability is transferred to two symmetric stable foci. 

At r = rhom ~ 13.9162 the unstable separatrices of the saddle return to the origin thus forming a homoclinic 
butterfly. This causes a " homoclinic explosion" in the phase space of the model that becomes filled in at 
once with countably many saddle periodic orbits that would further compose the skeleton of the Lorenz 
attractor. 

For rhom < r < rhet ~ 24.0579 the model exhibits transient chaos with two ultimate attractors: stable foci 
Oi 5 2- Such transient chaos is associated with a pre-turbulence regime. 

The value r = Th e t corresponds to the onset of the Lorenz attractor coexisting with stable foci 0\ and (92- 
The attraction basins of are bounded by the 2D cylindrically-shaped stable manifolds of two saddle 
periodic orbits that have earlier bifurcated from the homoclinic loops of the saddle at rh om ~ 13.9162. 
As r increases to tah ~ 24.7368 the saddle periodic orbits shrink the attraction basins and collapse onto 
the stable foci through a subcritical Andronov-Hopf bifurcation. 

For tah < t < ?~T ~ 31 the Lorenz equation possesses a genuinely strange chaotic attractor, known as the 
Lorenz attractor, containing no stable orbits. 

2.2. Canonical 2D bifurcation diagram of the Lorenz equation 

The pilot study of the dynamics of the Lorenz equation needs to be further enhanced by the bi-parametric 



examination of the model, including at large parameter values |Barrio fc Serrano , |2007 , |2009|. We w ill start 
off with the canonical bifurcation diagram, shown in the left panel of Fig. [3] (courtesy of |Shilnikov[ [T980 ) , 
of the Lorenz equation that depicts the basic bifurcation curves in the (r, a)-parameter plane with fixed 
6 = 8/3. The right panel of Fig. [3] sketches en route fragments of the formation of the Lorenz attractor 



on the pathway, a = 10 |Afraimovich et ah , 1977; Kaplan &; Yorke, 1979| through the bifurcation curves. 
For r < 1, Eq. ([!]) has a single stable equilibrium state at the origin. This equilibrium state undergoes 
a pitch- fork bifurcation at r = 1, so that for r > 1 the origin becomes a saddle, O, of the topological 
type (2,1) due to the characteristic exponents A3 < A2 < < Ai. The saddle has a ID unstable manifold, 
Wq is made of O itself and a pair of ID unstable separatrices, Ti and T2 (due to Ai) entering the saddle 
as t ^ — 00, and a 2D stable manifold, Wq, containing the leading (due to A2) invariant z-axis; the 
eigenvector due to A3 determines the non- leading or strongly stable direction, Wq in Wq. After the pitch- 
fork bifurcation, the separatrices of the saddle tend to two symmetric attractors - equilibrium states, 



01,2 — V — ±\/6(r — 1), z = r — 1) (Fig. [3j^A)) that become the global attractor for all trajectories in 
the phase space of the Lorenz equation other than in Wq. 

A homoclinic butterfly bifurcation occurs in the Lorenz equation when both separatrices, Ti and T2, 
of the saddle come back to the origin along the z-axis (Fig. |3^B)). In virtue of the symmetry of the 
Lorenz equation, such homoclinic loops are always formed in pairs, and therefore constitute bifurcations 
of codimension-one, in general. The bifurcation of the homoclinic butterfly takes place on the curve l\ in 
the (r, a)-parameter plane. Bifurcations of the separatrices of the saddle at the origin are crucial for the 
Lorenz attractor. 

The very first homoclinic butterfly made of two separatrices looping a single round about the equilib- 
rium states Oi ? 2, causes the homoclinic explosion in the phase space of the Lorenz equation. This bifurcation 
gives rise to a onset of countably many saddle periodic orbits that forming an unstable chaotic (saddle) set 
, which is not an attractor yet. For this explosion to happen, the so-called saddle value S — \\ + A2, which 
is the sum of the leading characteristic exponents of the saddle, must be positive; alternatively, the saddle 
index v = | A2I/A1 < 1. Otherwise, if S < [y > 1) the homoclinic butterfly produces a symmetric figure-8 
periodic orbit in the aftermath of the gluing bifurcation through which two stable periodic orbits merge 
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after flowing into the separatrix loops. L. Shilnikov Shilnikov, 1968| pointed out two more conditions, in 
addition to the primary one (1): a ^ 0, or v ^ 1, needed for the separatrix loop of the saddle in R 3 and 
higher dimensions to produce a single saddle periodic orbit; they are: (2) T E W ss , i.e. the separatrix 
comes back to the saddle along the leading direction, (3) the so-called separatrix value (A in the mapping 
Eq. Q below) does not vanish, its sign determines whether the separatrix loop is oriented or twisted, 
and hence the stable manifolds of the saddle periodic orbit are homeomorphic to a cylinder or a Mobius 
band. Otherwise, he predicted |Shilnikov , 1981] that the Lorenz attractor could be born right near such 
codimension-2 bifurcations, termed resonant saddle, orbit- and inclination- flip, correspon dingly |Shilnikov 



et al, 1998,2001 , as it occurs in the Shimizu-Morioka and similar models Shilnikov, 1986 Robinson. 1989 



Rychlic, 1990] below. 



Out of many saddle periodic orbits, which exploded the phase space of the Lorenz equation, two are the 
special, Li,2, ones as they demarcate the thresholds of the "interior" of the chaotic unstable set (Fig.[3jC)). 
After the homoclinic butterfly bifurcation, in the region between the bifurcation curves l\ and I2 in the 
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Figure 3. Left panel: the (r, cr)-parameter plane depicting the primary bifurcation curves and the stages of the formation of 
the Lorenz attractor that are sketched in the right panel. Curve l\ corresponding to the primary homoclinic butterfly shown 
in (B) ; I2 being the first boundary of the existence region of the Lorenz attractor: stages (C)-(D)-(E); £3 corresponding to a 
subcritical Andronov-Hopf bifurcation; and 1 4 and Z5 (F), corresponding to the homoclinic loops (depicted in in the insets) 
with symbolic kneadings {+1, —1,0} and {+1, — 1, +1, 0}, respectively. Courtesy of [Shilnikov} |1980| . 



(r, a)-parameter plane, the separatrices Ti and Y2 of the saddle switch the targets: now the right/left 
separatrix tends the left/right stable focus 02,1- 

In order for the unstable chaotic set to become the Lorenz attractor, it must become invariant, i.e. a 
closed set containing all a-limit orbits, and hence no loose of trajectories escaping to stable foci Oi ; 2- This 
occurs on the bifurcation curve, I2, in the parameter space (Fig. [3^D)). To the right of the curve, the basin 
of the Lorenz attractor is shielded away from those of the stable equilibrium states by the 2D cylinder- 
shaped stable manifolds of the two "threshold" saddle orbits, that have simultaneously emerged from 
both separatrix loops, Ti 5 2 at the homoclinic explosion on the curve l\ in the parameter plane. 

As one moves further to the right in the parameter plane, the saddle orbits, Li 5 2 5 keep narrowing the 
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Figure 4. 3D version of Inset [3p). The birth of the Lorenz attractor (grey): the attraction basins of the stable foci (purple 
dots) being blocked away from the extreme separatrices (blue orbits), l\2, of the saddle, O, at the origin and other trajectories 
on the Lorenz attractor by the cylinder-shaped 2D stable manifolds Wf J1 (dark blue) of the saddle periodic orbits 1,2,1 in 
the phase space of the Lorenz equation. 



attraction basins of the foci Oi ? 2, and on the bifurcation curve 1% they collapse into the stable equilibria. The 
equilibrium states become saddle- foci of the (l,2)-type through a subcritical Andronov-Hopf bifurcation 
|Roshchin, 1978|. The topological (l,2)-type means that each saddle- focus has 2D unstable and ID stable 



manifolds; the latter is formed by two incoming separatrices. Some local properties of the saddle-foci can 
be revealed without evaluating their characteristic exponents explicitly. Let Ai < stand for the real 
stable exponent of Oi 5 2, and A2,3 stand for a comp_lex conjugate pair such that ReA2,3 > 0. Observe that 
the divergence of the vector field defined by Eqs. 
Ai + 2ReA 2 ,3 < 0. This implies Ai + ReA 2 ,3 < 0, i.e. 



is given by [—a — 1 — 8/3], which equals J2i=i ^ = 
the complex conjugate pair is closer to the imaginary 



axis than the real negative exponent, and hence the saddle- foci meet the Shilnikov condition [Shilnikov 



1965, 1967; Shilnikov & Shilnikov, 20071. Therefore, as soon as the saddle- focus possesses a homoclinic 



loop, such a bifurcation causes the abundance of periodic orbits nearby. Those periodic orbits constantly 
undergo saddle-node and period doubling bifurcations as the parameters are varied. Moreover, since the 
divergence of the vector field of the Lorenz equation is always negative, saddle-node bifurcations give rise 
to stable periodic orbits near the homoclinic saddle-focus bifurcation. Under fulfillment of some global 
conditions, a single Shilnikov saddle- focus bifurcation can lead to the formation of a spiral or screw-like 
attractor. However, a strange attractor due to the Shilnikov saddle-focus in a 3D system with a negative 
divergence is no genuinely chaotic set in the sense that it contains stable periodic orbits within. Because 
of that, such a chaotic attractor is called quasi- attractor thus referring to that because besides stable 
periodic orbits with weak basins it may have structurally unstable or non-transverse homoclinic orbits 



Afraimovich & Shilnikov, 1983 Shilnikov 



1994 



1997 . Note that systems in higher dimensions can possess 



genuinely strange attractors with the Shilnikov loop without stable periodic orbits, the so-called wild 
chaotic attractors |Shilnikov[ |1994[ |1997[ |Turaev fc Shilnikov[ |1998[ |Shilnikov[ |2002[ . 

The Lorenz attractor is non-hyperbolic because it includes the singularity at the origin - the saddle 
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Figure 5. (A) Heteroclinic connection (in dark color) between the saddle at the origin and two saddle- foci (blue spheres) 
being overlaid with the Lorenz attractor (green light color) on the background at the primary T-point (r = 30.38, a = 10.2). 
Orange spheres on the butterfly wings indicating the turning points around the right and left saddle-foci define the kneading 
sequence entries, {±1}, respectively. (B) A typical time evolution of the x-coordinate of the right separatrix of the saddle. 



equilibrium state with an ID unstable manifold while all other saddle periodic orbits on the attractor 
have 2D stable and unstable manifolds. Moreover, the manifolds of those orbits in the Lorenz attractor 
(self) cross transversally in the 3D phase space thereby producing only structurally stable (transverse) 
homoclinic and heteroclinic trajectories. This condition imperative for the Lorenz attractor will not always 
hold for larger parameter values and lead to homoclinic tangencies of the manifolds which are followed by 
saddle-node bifurcations in the Newhouse regions |Gonchenko et al } 1993, 1996 of the parameter plane 
of the system. Thus, in order for the Lorenz attractor to be strange and chaotic with no stable orbits, it 
must not include the [homoclinic] saddle- foci, Oi 5 2, as well as contain only structurally stable homoclinic 
orbits due to transverse intersections of the manifolds of saddle periodic orbits. 



3. Kneading invariants 

Chaos can be quantified by several means. One customary way is through the evaluation of the topological 
entropy. The greater the value of topological entropy, the more developed and unpredictable the chaotic 
dynamics become. Another practical approach for measuring chaos in simulations capitalizes on evaluations 
of the largest (positive) Lyapunov exponent of a long yet finite-time transient on the chaotic attractor. 

After the stable foci have lost the stability through the subcritical Andronov-Hopf bifurcation, the 
Lorenz equation exhibits the strange attractor of the iconic butterfly shape. The wings of the butterfly are 
marked with two symmetric eyes containing the saddle-foci isolated from the trajectories of the Lorenz 
attractor. This attractor is structurally unstable |Guckenheimer fc Williams , 1979; Afraimovich et al , 1983| 



as it undergoes bifurcations constantly as the parameters of the Lorenz equation are varied. The primary 
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cause of structural and dynamical instability of chaos in the Lorenz equation and similar models is the 
singularity at the origin - the saddle with two one-dimensional outgoing separatrices. Both separatrices fill 
in densely two spatially symmetric, [(x, y, z) <-> (— x, —y. z)), wings of the Lorenz attractor in the 3D phase 
space (see Fig. [5]). The Lorenz attractor undergoes a homoclinic bifurcation when the separatrices of the 
saddle change the flip-flop pattern of switching between the butterfly wings centered around the saddle- 
foci. At such a change, the separatrices comes back to the saddle thereby causing additional homoclinic 
explosions in phase space |Afraimovich et al~ 1977; Kaplan &; Yorke , [1979] . 

The time progression of the "right" (or symmetrical "left") separatrix of the origin can be described 
geometrically and categorized in terms of the number of flip-flops around the equilibrium states 0\ and 
O2 in the 3D phase space of the Lorenz equation (Fig. [5]). Or, alternatively, can be reduced to the time- 
evolution of the x-coordinate of the separatrix, as shown in panel B of Fig. [5| The sign-alternation of 
the x-coordinate suggests the introduction of a {±l}-based alphabet to be employed for the symbolic 
description of the separatrix. Namely, whenever the right separatrix turns around 0\ or O2, we write down 
+1 or —1, respectively. For example, the time series shown in panel B generates the following kneading 
sequence starting with {+1, —1, —1, —1, +1, —1, —1, +1, —1, . . .} etc. 




-0.06-0.04-0.02 0.02 0.04 0.06 




-0.06 -0.04 -0.02 0.02 0.04 0.06 



Figure 6. (A) ID Lorenz mapping (geometric model) with the discontinuity corresponding to the saddle at the origin in 
the 3D phase space of the Lorenz equation. Shown are the forward iterates of the "right" separatrix, + , of the discontinuity 
point. Iterates on the right, x > 0, and left, x < 0, branches of the mapping are assigned to kneading symbols of +1, and —1, 
respectively. Here v — 0.65, A — 0.7 and \i — —0.06 in Eq. |2]). (B) Alterative cusp-shaped mapping as it would be for the 
z- variable of the separatrix at the turning points, given by ZmaxW — 5 on the Lorenz attractor. 



In what follows we will introduce and demonstrate a new computational toolkit for the analysis of 
chaos in the Lorenz-like models. The toolkit is inspired by the idea of kneading invariants introduced in 
Milnor & Thurston, 1988|. A kneading invariant is a quantity that is intended to uniquely describe the 



complex dynamics of the system that admit a symbolic description using two symbols, here +1 and — 1. 
The kneading invariant is supposed to depend monotonically on the governing parameter so that any two 
systems can be compared and differentiated, or equivalently, ordered in terms of > and <, by the kneadings. 
Two systems with the Lorenz attractors are topologically conjugate when they share the same kneading 
invariant, see [Guckenheimer fc Williams , 1979; Rand, 1978 Malkin, 1991; Glendinning &; Sparrow, 1993 



Tresser & Williams, 1993| and the references therein. Moreover, for the Lorenz-like systems, the topological 



entropy can be evaluated though the kneading invariants by reducing consideration to piecewise monotone 
mappings of a closed interval |Glendinning fc Hall, 1996; Li &; Malkin, 2003 . 



Such mappings are closely related to the geometric models of the Lorenz attractor which are ID and 2D 
Poincare return mappings defined on some cross-section transverse to trajectories of the Lorenz attractor. 
The basic idea behind either geometric model capitalizes on extended properties of the local Poincare 
mapping near a homoclinic butterfly of the saddle |Guckenheimer fc Williams , 1979; Afraimovich et al. 
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1983; Shilnikov et a/., 1998,2001] . The mapping is assumed to meet a few other conditions related to the 
global properties of the flow far from the homoclinic butterfly. Such a ID constrained mapping shown in 
Fig. [6^A) can be written as: 

T: en+i = (M + A|e„r)-sign(£ n ), (2) 

here v — | A2 1/ Ai < 1 is the saddle index, fi controls the distance between the returning separatrices, I^i, 
and the saddle, and A is a scalar whose sign determines whether the homoclinic loops at \± — are oriented 



if A > , or twisted when A < 0, see |Shilnikov et ah , 1998,2001] for more details. 

Loosely speaking, this geometric model (|2j) should have no critical point on bo th branches, and moreove r 



Afraimovich et al. 



1983 



possess a property of strong stretching with a rate of expansion more than y/2 
This would guarantee that the Lorenz attractor will be densely filled in by the forward iterates of the 
separatrices, ± with no holes - lacunas containing isolated periodic orbits inside, stable or not. Strong 
stretching would also imply a monotone dependence of the kneading invariants on a governing parameter. 

In a symmetric system with the Lorenz attractor, the kneading invariant is assigned to quantify 
the symbolic description of either separatrix; in the asymmetric case one should consider two kneading 
invariants for both separatrices of the saddle. Thus, in respect it reflects quantitatively a qualitative change 
in the separatrix behavior, such as flip-flopping patterns, as the parameter of the system is changed. 

By construction, kneading invariants are proposed to serve as moduli of the topological equivalence 
that are employed to compare or contrast between any two Lorenz attractors or, equivalently, any two 
Lorenz-like systems. Due to the symmetry of the Lorenz mapping = T(£ n ) = T n (^o) from Eq. ([2]), 
forward iterates of the right separatrix, + , of the discontinuity point (resp. the saddle) are detected to 
generate a kneading sequence {ft n (0 + )} defined by the following rule: 



+1, if T n {0+) > 0, 
,(0 + ) = { -1, if T n (0+) < 0, 
0, if T n (0+) = 0; 



(3) 



here T n (0 + ) is the n-th iterate of the right separatrix + of the origin. The condition T n (0 + ) = is 
interpreted as a homoclinic loop, i.e. the separatrix returns to the origin after n steps. 



{+V1r1rH,+1,+1,+1,...,+1,...} 




{+VVV1,-1,-VW1...} 



-20 -15 -10 



5 10 15 20 25 




Figure 7. Truncated kneading sequences generated by the right outgoing separatrix of the saddle at the origin in the Lorenz 
equation at two distinct values of the parameters. 



The kneading invariant for the separatrix is defined in the form of a formal power series: 

00 

P(q) = Y J K n q n . (4) 



=0 



Setting q E (0, 1) make the series Q convergent. The smallest zero, if any, of the graph of Q in the 
interval q E (0, 1) yields the topological entropy, h(T) = ln(l/g*), of the ID mapping Q. 



10 Barrio, Shilnikov and Shilnikov 



Let us next draw a parallelism between the geometric model and the Lorenz equation: the kneading 
sequence {k u } comprised of only +ls of the mapping ^ corresponds to the right separatrix converging to 
the stable equilibrium state, 0\ (or possibly, a periodic orbit with x(t) > 0). The corresponding kneading 
invariant is maximized at {P m ax(q)} = 1/(1 — q)- When the right separatrix converges to an cj-limit set 
with x(t) < 0, like the left stable focus, O2 then the kneading invariant is given by {P m in(^)} = 1 — q/0- — q) 
because the first entry +1 in the kneading sequence is followed by infinite —Is. Thus, [{P m in(^)} 5 {-Fmax (<?)}] 
yield the range of the kneading invariant values; for instance [{P m in(l/2)} = 0, {P m ax(l/2)} = 2]. Two 
samples of the separatrix pathways shown in Fig. [7] generating the following kneading invariants 

P A (l/2) = +1 - 1/2 - 1/4 + 1/8 + 1/16 + 1/32 + 1/64 . . . + l/2 n ... = 1/2, 
P B (l/2) = +1 - 1/2 - 1/4 - 1/8 - 1/16 - 1/32 - 1/64 ... - l/2 n . . . = 0, 

illustrate the parallelism. 

To conclude this section we point out that there is another approach for constructing ID return 
mappings through the evolution of the z-variable of the separatrix of the Lorenz equation. The mapping 
generated by the turning points where 4 ax (^) = on the attractor has a distinct cusp-shaped form depicted 
in Fig. [6](B) . The point corresponding to the cusp is used for the initiation of the kneading sequence. Due 
to its unimodal graph with two, increasing and decreasing, segments the corresponding formal power series 
is then defined as follows: 

00 n 

p (q) = ^2K n q n , where K n = JJft*, (5) 

71=0 2=0 

i.e. Hi n — Hi n • K n —\. 



4. Kneading scanning of the Lorenz equation 

In this section we carry the concept of the kneading invariants over to numerical studies of fine structures of 
chaos in the Lorenz equation. For sake of simplicity in this pilot phase, we employ a rough idea for defining 
the kneading sequences of +ls and —Is that relies on whether the right separatrix of the saddle makes a 
revolution around the right equilibrium state, Oi, or the left one, O2, respectively, in the (x, z)-projection of 
the 3D phase space. One can utilize even a simpler approach where the sign of the current kneading entry 
in the sequence is determined by the sign of the x-coordinate of the separatrix at the off-lying turning 
points, max|x(t*)| (on the butterfly wings), given by x'(t*) = 0, see the trace in Fig. [5^B). There are 
other alternative ways for defining kneading entries, involving cross-sections, z\t) etc, that are not free of 
certain limitations either. In the future we plan to enhance the current kneading algorithms by utilizing 
the incoming ID separatrices of the saddle-foci and finding the winding numbers around them instead. 

In this computational study of the Lorenz equation and two other models below, we consider partial 
kneading power series truncated after the first 50 entries: P$o(q) = Yln°=o ^nq 71 - The choice of the number 
of entries is not motivated by numerical precision, but by simplicity, as well as by a resolution of the bitmap 
mappings for the bi-parametric scans of the models. One has also to figure the proper value of q: setting it 
too small makes the convergence fast so that the tail of the series has a little significance and hence does 
not differentiate the fine dynamics of the Lorenz equation on longer time scales. 

At the first stage of the routine, we perform a bi-parametric scan of the Lorenz equation within a specific 
range in the (r, cr)-plane. The resolution of scans is set by using mesh grids of [1000 x 1000] equa lly- dist anced 
points. Next by accurately integrating the separatrix using Taylor series software TIDESQ [ Abad et al 



2011a|b ; Barrio et al } 2011c we identify and record the sequences {/%}5o for each point of the grid in 



the parameter plane. Then we define the bi-parametric mapping: (r, a) —> Pso(q) with some appropriately 
chosen g, the value that determines the depth of the scan. The mapping is then colorized in Matlab by 
using various built-in non-linear spectra ranging between to P^ in and P^f 1 *, respectively. In the mapping, 
a particular color in the spectrum is associated with a persistent value of the kneading invariant on a level 
curve. Such level curves densely foliate the bi-parametric scans. 



freeware TIDES is a versatile numerical ODE solver for integration of ODEs with an arbitrary precision, especially for chaotic 
systems. 
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Figure 8. (A) Kneading-based bi-parametric scan revealing multiple T-points and saddles that organize globally complex 
chaotic dynamics of the Lorenz equation in the (r, a) parameter plane. Solid-color regions associated with constant values of 
the kneading invariant correspond to simple dynamics dominated by stable equilibria or stable periodic orbits. The borderline 
between the brown and blue region corresponds to the bifurcation curve of the homoclinic butterfly. The border line between 
the blue and yellow-reddish region corresponds to the formation of the Lorenz attractor (below a c± 50). (B) Zoom of the 
vicinity of the primary T-point at (r = 30.4, a = 10.2) to which a homoclinic bifurcation curve spirals onto. Data for the 
homoclinic curves (in blue) are courtesy of Yu. Kuznetsov. (C) Original bifurcation diagram of the Lorenz equation depicting 
the two detected T-points and primary homoclinic bifurcation curves; courtesy of [Bykov| [1993 
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Figure 9. The (r, a)-bifurcation diagram of the Lorenz equation depicting the existence region (shaded) of the Lorenz 
attractor. The second bifurcation curve, passing through the primary T-point, Q(r = 30.4, a = 10.2), crosses the first 
boundary, l a (h hi Fig.|3|, at a point, labeled by K, thus closing the existence region. Crossing the branch l a rightward above 
the point K does lead to the emergence of the Lorenz attractor: after a long chaotic transient, the separatrix of the saddle 
will be attracted to either stable foci, 0\ or O2. Courtesy of [Bykov &; Shilnikov 1992 . 



Figure [8] represents the kneading-based color scan of the dynamics of the Lorenz equation mapped 
onto a fragment of the (r, a)-parameter plane. In the scan, a window of a solid color corresponds to a 
constant kneading invariant. In such windows the dynamics of the Lorenz equation are dominated by 
simple attractors such as stable equilibria or stable periodic orbits to which the separatrix converge. A 
quick examination of the kneading definition Q reveals that the kneading invariant does not vary at a 
supercritical Andronov-Hopf bifurcation and a pitch-fork bifurcation describing continuous transitioning 
between stable symmetric and asymmetric periodic orbits. This holds true for a period-doubling bifur- 
cation too, because the kneading sequence, say {(+1, — 1, — I) 00 }, inherits the same block repeated twice 
{(+1, —1, —1, +1, —1, — I) 00 } in the code for the periodic orbit of a double period and so forth. While the 
kneading technique does not detect such safe bifurcation boundaries |Shilnikov et aZ. , 1998,2001 having 
crossed through which the phase point does not run far from an old attractor to a new one, though it 
detects dangerous boundaries well, including a generic saddle-node bifurcation, homoclinic bifurcations, 
and some other. 

A borderline between two solid-color regions corresponds to a bifurcation through a dangerous bound- 
ary which is associated with a jump between the values of the kneading invariant. For example, the 
borderline in Figure [8] between the brown region with the kneading sequence {(+I) 00 } and the blue region, 
with the kneading sequence {+1, ( — I) 00 }, corresponds to the primary homoclinic butterfly of the saddle. 
The second borderline of the blue region corresponds to the onset of the Lorenz attractor existing on the 
right from it. One can see that above a ~ 50 this border is adjoined by windows of solid colors thus 
indicating that the separatrices start converging to stable equilibria after chaotic transients, long or short. 
Indeed, crossing the curve, I2 (or l a in Fig. [3]), above a = 18 does not imply that the Lorenz equation will 
have a strange attractor, but a chaotic saddle set |Bykov fc Shilnikov] |1989 , 1992| . 

What the proposed kneading technique does extraordinary well, compared to the bi-parametric screen- 
ing based on the finite-time Lyapunov exponent approach (see below), and what we developed it for, is 
the detection of bifurcations within the Lorenz strange attractor. The corresponding yellow-reddish regions 
in the parameter plane in Fig. [8] clearly demonstrate the evidence of the parametric chaos that, like in 
turbulence, enriched by vortices of multiple T-points. Panel B of this figure depicts the kneading mapping 
near the left-bottom corner of the bifurcation diagram in panel A. In it, the black (blue in panel A) region 
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corresponds to the chaotic saddle dynamics transitioning into the Lorenz attractor after the mapping color 
shifts to the yellow-reddish spectrum. Blue curves in panel B are the bifurcation curves of some homoclinic 
orbits with short admissible kneadings. One can see from this panel that the mapping spectrum is clearly 
foliated by the kneading invariant level curves of the colors gradually progressing from red to yellow. This 
indicates that the new born Lorenz attractor, while being structurally unstable and sensitive to parameter 
variations, persists initially the pseudo-hyperbolic property because the foliation remains uniform, and 
transverse to the classical pathway a — 10 (Fig. [3]). The homogeneous foliation starts breaking around a 
saddle point after which one singled bifurcation curve spirals onto the primary T-point. Far from this point, 
the curve corresponds to the formation of the homoclinic loop with the kneading (1,-1,-1,-1,0); i.e. the 
right separatrix makes one excursion around the saddle- focus Oi, followed by three revolutions around 
the saddle- focus O2, and then returns to the saddle at the origin. While moving along the spiraling curve 
toward the T-points, the separatrix makes progressively more turns around O2, or more precisely around 
the ID incoming separatrix of this saddle- focus. With each incremental turn around O2, the separatrix 
comes closer to O2 while the bifurcation curve becomes one scroll closer to the T-point simultaneously. 
Due to this feature the T-point Qo(r = 30.4, a — 10.2) is called a Terminal point. The T-point corresponds 
to the following symbolic sequence {+1, (— I) 00 }. In virtue of the symmetry of the Lorenz equation, the 
T-point actually corresponds to a closed heteroclinic connection involving all three saddle-equilibria, see 
Figs. [I] and [5j The merger of the right (left) separatrix of the saddle O with the incoming separatrix of the 
saddle- focus O2 (Oi), increases the codimension (degeneracy) of this heteroclinic bifurcation to two; note 
that intersections of the 2D unstable manifolds of the saddle-foci, with the 2D stable manifold of the saddle 
at the origin are transverse in the 3D phase space in general. Breaking the ID heteroclinic connection gives 
rise to a primary homoclinic orbit to the saddle-focus, as well as to a heteroclinic connection between 
both saddle- foci (see the sketch of the bifurcation unfolding of the T-point in Fig. [I]). The corresponding 
bifurcation curves of these homoclinic and heteroclinic bifurcations originating rightward from the T-point 
bound a sector containing subsidiary T-points |Bykov, 1993 Glendinning fc Sparrow[ [l984| . Each new T- 
point produces other self-similar structures scaled like fractals. Panel C shows two such identified T-points: 
primary Qo and secondary Q\ (r = 85. a = 11.9). The primary codimension-two T-point in the Lorenz 
equation was originally discovered by Yudovich |Pertovksaya fc Yudovich , 1980] . 
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Figure 10. Finite-time largest-Lyapunov exponent, Lmax, scan of the Lorenz equation showing no sign of spiral structures 
in the (r, cr)-parameter plane. The dark region corresponds to trivial attractors, where L max < 0, while the red color indicates 
L max > in chaotic regions. The red dot points out the location of the primary T-point. 
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As soon as the saddle-foci and their bifurcations become involved in the dynamics of the Lorenz 
equation near the primary T-point, the Lorenz attractor loses the purity of the genuine chaotic attractor 
that used to have neither stable periodic orbits nor non-transverse homoclinic trajectories. It transforms 
into a quasi-chaotic attractor with weakly stable orbits of small basins and nontransverse homoclinic orbits. 
The idea of non-transversality or tangency was employed in |Bykov fc Shilnikov , 1989, 1992| to numerically 
identify the second boundary, Ik , in addition to the first one, l a (J2 in Fig. |3j), that bounds the existence 
region of the Lorenz attractor in the parameter plane, see Fig. |9j Note that Ik crosses the initial boundary, 
l a . This means that above the intersection point, crossing l a {I2 in Fig. [3]) rightwards does not guarantee 
that the basin of the Lorenz attractor will necessarily be isolated from the basins of the stable foci, 
by the stable manifolds of the saddle periodic orbits, (Fig. [4]). This implies that the separatrices of 
the saddle will demonstrate chaotic transient behavior, long or short, prior to them converging to Oi 3 2, 
Shilnikov] [19861 [l989| |Bykov fc Shilnikov[ [1989} |Shilnikov] [1991] |Bykov fc Shilnikov] [l992j |Shilnikov 



see 



1993; Shilnikov et aL 1993 for full details. 



The other feature of the boundary, Ik, is that it passes through the primary T-point, thereby separating 
the existence region of the Lorenz attractor from bifurcations of the saddle- foci, and consequently from all 
subsidiary T-points existing on the right from it in the parameter plane. Here the chaotic dynamics of the 
Lorenz equation become even wilder and less predictable |Shilnikov , 1994, 1997; Turaev & Shilnikov, 1998 



Shilnikov, 2002| . The indication to that is panel A of Fig. |8j that reveals, through the kneading scan, a 



parametric turbulence in the (r, cr)-parameter plane with fractal explosions in the forms of multiple spiral 
structures - "tornado eyes" centered around T-points. Note that basins of spiraling T-points are separated 
by corresponding saddles. One can spot self-similar smaller-scale spiral structures within large-scale ones 
and so forth. The richness of such fractal structures in the parameter plane results from the synergy of the 
Lorenz-like dynamics amplified by chaos induced by the Shilnikov saddle-foci. 

To conclude this section we contrast the scans of the Lorenz equation obtained using the proposed 
kneading technique with the sweeps based on the evaluation of the largest Lyapunov exponent, L max , for 
the separatrices of the saddle evaluated over a finite time interval |Barrio fc Serrano} |2007[ 2009; Barrio 



et aL 2011a 



Figure 10 shows a fragment of the typical bi-parametric sweep of the Lorenz equation: 
the dark region at small values of the Reynolds number, r, is where L max is negative on the stable foci. 
The sweep yields a clear borderline between the regions of the simple and chaotic attract ors. The chaotic 
region (yellow-reddish) is characterized by a small positive Lyapunov exponent, variations of which are 
not significant enough to reveal fine structures, as spiraling T-points. The method can detect well stability 
islands in the biparametric diagram which correspond to emergent stable periodic orbits. 

Two panels in Fig. 11 represent the expanded scans of the of Lorenz equation: panel (A) is the kneading 
invariant mapping on the grid of [1000 x 1000] points and panel (B) shows L maar based sweeping of the 
(r, a)-parameter plane. While panel A reveals a chain of large-scale spiral hubs around T-points, panel B 
reveals none but stability windows (dark). We would like to point out that the stability windows can also 
be detected in the kneading scan in panel A showing the border of the solid color (dark yellow) island 
stretched horizontally at small values of the parameter a. 



5. The Shimizu-Morioka model 

Next we perform the kneading-based bi-parametric scanning of another classical three-dimensional system 



called the Shimizu-Morioka model [Shimizu fc Morioka[ |1980[ |Shilnikov[ |1986[ |1989[ |1991[ |1993| : 



x 



y, y 



Ay 



xz, 



-OLZ + X 



(6) 



here, a and f3 are positive bifurcation parameters. Like the Lorenz equation, this Z2-symmetric model has 
3 equilibrium states: a saddle of the (2,l)-topological type at the origin, and two symmetric stable-foci or 
saddle- foci of the (l,2)-topological type. 

This model was originally introduced to examine a pitch-fork bifurcation of the stable figure-8 periodic 
orbit that gives rise to multiple cascades of period doubling bifurcations in the Lorenz equation at large 
values of the Reynolds number r. It was proved in | Shilnikov et al. , 1993] that the Eqs. ^ are a universal 
normal form for several codimension-three bifurcations of equilibria and periodic orbits on Z2-central man- 
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Figure 11. (A) Kneading-based scan revealing a fractal structure and a chain of spiral vortices centered at T-points alternating 
with saddles in the extended (r, a) -region of the Lorenz equation. (B) Lmax - based sweeping of the Lorenz equation singes out 
stability windows (dark) corresponding to steady state and emergent periodic attractors with L max < within the chaotic 
region (white) associated with L m ax > 0. 




Figure 12. A partial (a, A)-diagram of the Shimizu-Morioka model depicting initial bifurcation curves and corresponding 
insets for the separatrix behaviors. Legend: AH stands for a supercritical Andronov-Hopf bifurcation, HB stands for the 
homoclinic butterfly made of two separatrix loops; SN stands for a saddle-node bifurcation of periodic orbits; it connects 
the codimension-two points, the resonant saddle a — on HB and the Bautin bifurcation at GH. LA stands for the Lorenz 
attractor formation; A — stands for an orbit-flip bifurcation for the double-loop homoclinics on H2. The dashed line separates, 
with good precision, the Lorenz attractor region from the region of a quasi-attractor (below). Vertical Pathway showing the 
gluing bifurcation on HB. Courtesy of iShilnikov et al. 1998,2001 . 



ifolds. The model turned out to be very rich dynamically: it exhibits various interesting global bifurcations 
|Shilnikov , 1986, 1991, 1993| including T-points for heteroclinic connections. 
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While the model inherits all basic properties of the Lorenz equation, in addition, and of special interest, 
are two homoclinic bifurcations of codimension-two: resonant saddle with the zero saddle value or the saddle 
index v = 1, and the orbit-flip bifurcation where the separatrix value A in Eq. ^ vanishes. Recall that 
the sign of the separatrix value determines whether the homoclinic loop, here double-pulsed, is oriented 
or flipped. These codimension-two points globally organize the structure of the compact (a, A)-parameter 
region of the Shimizu-Morioka model, including structural transformations of the Lorenz attractor in the 
model, including the emergence of stability islands - lacunae inside the attractor. 

Figure 13 represents a partial (a, A)-diagram of the Shimizu-Morioka model and depicts primary 
bifurcation curves and the corresponding phase portraits of the separatrix behaviors. Among those are the 
Andronov-Hopf bifurcation curve, AH, above which the equilibrium states, are stable, and saddle- 
foci below. This bifurcation is primarily supercritical, but becomes subcritical at smaller values of the 
parameter a. The bifurcation curve, HB, corresponds to the formation of the homoclinic butterfly, or 
flgure-8 at larger values of a. The transition between the branches of the curve is no bifurcation like the 
inclination-switch because the saddle value, a is negative here. The pathway at a = 1.15 demonstrates 
the evolution of the simple Morse-Smale dynamics of the model from the stable equilibrium states to a 
stable symmetric periodic orbit. This orbit emerges through a gluing bifurcation after that two stable 
periodic orbits existing below the supercritical Andronov-Hopf curve AH form the homoclinic butterfly 
with a < 0. The saddle value becomes positive to left from the codimension-two point, labeled by a = 
on the bifurcation curve HB. The left segment of the curve HB is similar to the bifurcation curve, Zi, of 
the Lorenz equation (Fig. [5]); namely, the homoclinic butterfly with a > causes the homoclinic explosion 
in the phase space of the Shimizu-Morioka model. The curve LA [LA\ in Fig. 13) being an analog of the 
curve I2 in the bifurcation diagram in Fig. [3] is the upper boundary of the existence region of the Lorenz 
attractor is the parameter plane of this model. Below LA the separatrices of the saddle no longer converge 
to stable periodic orbits but fill in the strange attractor, like in Fig. |1J The bifurcation unfolding of the 
homoclinic resonant saddle includes various bifurcation curves: among them in the figure we depict, in 
addition to LA, the curve, SN, corresponding to saddle-node of merging stable (through the supercritical 
Andronov-Hopf bifurcation) and saddle (through the homoclinic bifurcation on HB) periodic orbits. The 
curve labeled by H2 corresponds to a pair of the double-pulsed homoclinic loops. Continued further away 
from the codimension-two point, a = 0, the curve H2 frames the chaotic region in the parameter plane. 
The point, A = 0, on this curve is a codimension-two orbit-flip homoclinic bifurcation: to the left of it, the 
loops become flipped like the median of a Mobius band. The dashed line, AZ originating from the point 
A = is alike the curve in Fig. [9] for the Lorenz equation. This curve is the second boundary of the 
existence region (above) of the Lorenz attractor in the Shimizu-Morioka model, that separates wild chaos 
with homoclinic tangencies (below). This like passes through the primary T-point in the (a, A)-parameter 
plane. Figure [13] depicts a few more bifurcation curves originating from the point, A = 0: two four-pulses 
homoclinic curve terminating at subsidiary T-points and the curve labeled by "—1" corresponding to a 
period doubling bifurcation. An intersection of the dashed line with a homoclinic bifurcation curve (see 
Fig. 13) corresponds to another orbit-flip bifurcation and so forth. 

Indeed the skeleton of the bifurcation set of the Shimizu-Morioka is more complex. The detailed bifur- 
cation diagram is shown in the top panel of Fig. 13 It reveals several T-points, the pitch-fork bifurcation 
curve, PH , among other bifurcation curves for various homoclinic and heteroclinic connections. The de- 
tailed description of the bifurcation structure of the Shimizu-Morioka model is out of scope of this paper. 
The curious reader can find a wealth of information on bifurcations of the Lorenz attractor in the original 
papers Shilnikov] 1 1986 , 1989, 1993; Shilnikov et al. , 1993| . We point out that those bifurcation curves were 
continued in the (a, A)-parameter plane of the model using A. Shilnikov's home-made software also based 
on the symbolic kneading toolbox. 

The bottom panel of Fig. 13 is a de-facto proof of the new kneading invariant mapping technique. The 
panel represents the bi-parametric scan in color of the dynamics of the Shimizu-Morioka model that is 
based on the evaluation of the first 50 kneadings of the separatrix of the saddle on the grid of 1000 x 1000 
points in the (a, A)-parameter region. Getting the mapping took about one hour on a high-end workstation 
without any parallelization efforts. The color scan reveals a plethora of large-scale T-points, as well as 
nearby smaller ones (Fig. 14) invisible in the given parameter range, as well as the saddles separating 
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Figure 13. Detailed (a, A)-parameter plane of the Shimizu-Morioka model obtained by the parameter continuation method 
(courtesy of Shilni kov et al\ |1993] ) and by the bi-parametric scan based on the kneading invariants. The scan revealing 
multiple T-points and saddles that globally organize complex chaotic dynamics of the model. Solid-color regions associated 
with constant values of the kneading invariant correspond to simple dynamics dominated by stable equilibria (brown) or stable 
periodic orbits (blue). The border between the brown and blue regions corresponds to the bifurcation curve of the homoclinic 
butterfly. The codimension two point, a = 0, gives rise to loci of bifurcation curves including LA\ below which the Lorenz 
attractor exists. Bifurcation loci of the other codimension two point, A — (orange zone) giving rise to subsidiary orbit-flip 
bifurcations on turns of spirals around T-points, are separated by saddles (two large scale ones) in the parameter plane. 
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Figure 14. Zoom of the (a, A)-parametric mapping in Fig. p"3|B) near the primary T-point revealing self-similar structures 
embedding smaller-scale spirals around secondary T-points in the Shimizu-Morioka model. 



spiral structures. 

The solid-color zones in the mapping correspond to simple Morse-Smale dynamics in the model. These 
trivial dynamics are due to either the separatrix converging to the stable focus 0\ (O2) and emergent 
periodic orbit with the same kneading invariant (brown region), or to the symmetric and asymmetric stable 
figure-8 periodic orbits (dark blue region). The borderlines between the simple and complex dynamics in 
the Shimizu-Morioka model are clearly demarcated. On the top it is the curve, LAi, (see the top panel 
of Fig. 13). The transition from the stable 8-shaped periodic orbits to the Lorenz attractor (through the 
boundary, LA2) is similar though more complicated as it involves a pitch- fork bifurcation and bifurcations 
of double-pulsed homoclinics, see |Shilnikovj |1993[ Shilnikov et al\ |1993| for details. 

One can clearly see the evident resemblance between both diagrams found using the bifurcationaly 
exact numerical methods and by scanning the dynamics of the model using the proposed kneading invariant 
technique. The latter reveals a richer structure providing finer details. The structure can be enhanced 
further by examining longer tails of the kneading sequences. This allows for the detection of smaller-scale 



spiral structures within scrolls of the primary T- vortices, as predicted by the theory Bykov, 1993 . Figure 14 



shows a magnification of the scan of the Shimizu-Morioka model near the primary T-point that depicts 
several other small-scale T-points. 

Finally, we compare the new kneading scanning apparatus with the customary bi-parametric sweeping 
(shown in Fig. 15) of the Shimizu-Morioka model that is based on the evaluation of the Lyapunov exponent 
spectrum computed over a finite time interval |Gonchenko et aT , 2005| . Likewise the case of the Lorenz 
model, the sweeping based on the Lyapunov exponents shows no sign of spiral or saddle structures inside 
the region of deterministic chaos. The regions of the solid colors are associated with the sign of the largest 
Lyapunov exponent, L m ax- negative L max values correspond to steady state attractors in the green region; 
L max = corresponds to periodic attractors in the blue region; and L max > is associated with chaotic 
dynamics in the model in the region. Note blue islands in the red-colored region that correspond to stability 
windows in chaos-land. In such windows the Lorenz attractor has an emergent lacuna containing, initially, 
a single symmetric saddle periodic orbit. The orbit then undergoes a pitch- fork bifurcation that makes it 
stable. The basin of the stable orbit, which is first bounded by the 2D stable manifold of two asymmetric 
saddle periodic orbits, increases so that the stable orbit starts to dominate over chaotic dynamics in the 
corresponding stability window. These bifurcations underlie the transition from simple dynamics (blue 
region) due to the symmetric stable periodic orbit to chaos through the curve, H2, as the parameter a in 
decreased. 



Kneadings, symbolic dynamics and painting Lorenz chaos 19 




Figure 15. (A) Attractors of the Shimizu-Morioka model being differentiated by the sign of the largest Lyapunov exponent, 
Lmax- Color legend for the attractors of the model: green-stable equilibrium states, Lmax < 0; blue - stable periodic orbits 
with a nodal normal behavior, Lmax — 0; magenta - a periodic orbit with a focal normal behavior; red - chaotic attractors 
with Lmax > 0, with identified lacunae. Courtesy of [Gonchenko et aT\ |2005] . (B) Lorenz attractor (blue) with a lacuna 
containing a symmetric figure-8 periodic orbit (dark red). 



6. 6D optically pumped laser model 

The coexistence of multiple T-points and accompanying fractal structures in the parameter plane is a 
signature for systems with the Lorenz attractor. A question though remains whether the new computational 
technique will work for systems of dimensions higher than three. In fact, to apply the technique to a generic 
Lorenz-like system, only wave forms of a symmetric variable progressing in time, that consistently starts 
from the same initial condition near the saddle is required. Next is an example from nonlinear optics - a 6D 



model of the optically pumped, far infrared red three-level molecular laser jMoloney et a/. , 1989; Forysiak 
et al. , 1991] given by 



P 

P21 

P23 

P31 

D21 

D23 



-cr/3 + gp23, 

-P21 ~ PP31 + aD 2 l, 

-P23 + (3D 23 - ap3i? 

~P31 + PP21 + ap23, 

-b(D2i-D° 21 )-4ap2i-2/3 P 23, 
-b(D 2 3 - £>§ 3 ) - 2«P2i - 4/3^23- 



(7) 



Here, a and b are the Rabi flopping quantities representing the electric field amplitudes at pump and emis- 
sion frequencies. The parameter a is a natural bifurcation parameter as it is easily varied experimentally. 
The second bifurcation parameter, 6, can be varied to some degree at the laboratory by the addition of 
a buffer gas. This system presents, like the Lorenz equations, a symmetry (P,P2i,P23,P3i, ^2i 5 D23) <H> 
(— /3,£>2i 5 ~~ P23, — P31, D21, D23). The laser model has either a single central equilibrium state, O (with 
/3 = 0), or additionally, through a pitch- fork bifurcation, a pair of symmetric equilibrium states, Oi^2 (with 
P > 0); the stability of the equilibria depends on the parameter values. 

Optically pumped, far infrared lasers are known to demonstrate a variety of nonlinear dynamic behav- 



iors, including Lorenz-like chaos |Hubner et al. , 1995] . An acclaimed example of the modeling studies of 
chaos in nonlinear optics is the two level laser Haken model |Haken , 1975] to which the Lorenz equation 
can be reduced. A validity that three level laser models would inherently persist the Lorenz dynamics 
were widely questioned back then. It was comprehensively demonstrated in jForysiak et al. , 1991 in 1991 
that this plausible laser model possesses a plethora of dynamical and structural features of the Lorenz-like 
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Figure 16. (A) Lorenz attractor with a lacuna in the laser model at a = 1.14, b = 0.2, q = 50 and a = 10. (B) (a, &)- 
bifurcation diagram of the model for g — 52 and a = 1.5. BP and HB here denote the pitch-fork and Andronov-Hopf 
bifurcations, respectively. HO and HE denote the branches of the primary homoclinic (of the saddle) and heteroclinic orbits 
(of the saddle- foci). C2 is the codimension two Khorozov- Taken point for the equilibrium state with double zero eigenvalues, 
and T is the primary terminal point. The spiraling curve connects the T-point with the homoclinic resonant saddle on HO, 
near which separatrix loops are double pulsed ones. Courtesy of [Forysiak et al. 1991] . 



systems, including the Lorenz attractor per se (with lacunae as well), similar Andronov-Hopf, Z2 pitch- 
fork, various homoclinic and heteroclinic bifurcations including codimension-two ones T-points, see Fig. [TBI 
Similar structures were also discovered in another nonlinear optics model for a laser with a saturable ab- 
sorber which can be reduced to the Shimizu-Miorioka model near a steady state solution with triple zero 
exponents |Vladimirov fc Volkov , 1993 

The laser model (|7|) is quite rich in bifurcations; the list also includes a super-critical Andronov-Hopf 
bifurcation of the central equilibrium state that gives rise to a stable figure-8 shaped periodic orbit (in 
proper projections) for the parameter values to the left of the bifurcation curve, AHo, in the bifurcation 
diagrams shown in panels B and C of Fig. [TT| Observe from the diagram that the curve AHo originates 
from the point labeled, BT. This point corresponds to a codimension-two Khorozov- Takens bifurcation of 
an equilibrium state with two zero Lyapunov exponents. The bifurcation is an extension of the Bogdadov- 
Takens bifurcation on a symmetric central manifold. The unfolding of this bifurcation includes two more 
curves: AH\2 standing for the Andronov-Hopf bifurcation for the secondary equilibrium states, 0\^\ and 
Horn standing for a homoclinic connection made of two symmetric separatrix loops bi-asymptotic to 
the saddle, O. The continuation of the curve, Horn, away from BT reveals rather peculiar details that 
substantially organize the bifurcation diagram of this laser model. Near BT the curve, Horn, corresponds 
to a homoclinic figure-8 of the saddle with a negative saddle value, on the ^-symmetric 2D central 
manifold in the 6D phase space of the model. Recall that the figure-8 homoclinic connection stands for 
the case where the ID unstable separatrices come back to the saddle, O from the symmetrically opposite 
directions along the eigenvector corresponding to the leading stable exponent at the equilibrium state. 
This bifurcation gluing two stable periodic orbits emerging from through the supercritical Andronov- 
Hopf bifurcation gives rise to the stable symmetric figure-8 periodic orbit existing nearby BT. As the 
curve, Horn, is continued further away from BT, the stable leading direction at the saddle, O, changes: it 
becomes the invariant /3-axis (like the z-axis in the Shimizu-Morioka model) so that the separatrix loops 
start returning tangent to each other and hence forming the homoclinic butterfly. Nevertheless, this is a 
gluing bifurcation, not a codimension-two bifurcation of the change of the leading direction (inclination 
switch) as the saddle value remains negative on this branch of the curve, Horn. The saddle value vanishes, 
making the saddle resonant at the codimension-two point, and becomes positive further down on Horn. As 



Kneadings, symbolic dynamics and painting Lorenz chaos 21 




Figure 17. (A) Bi-parametric scan of the laser model featuring the T-points and saddles typical for the Lorenz-like systems, 
mapping the dynamics of the 6D optically pumped far-infrared-red laser model onto the (electric-field-amplitude, omission- 
frequency)-diagram at g = 50 and a — 1.5. Solid-color windows and fractal regions correspond to trivial and chaotic dynamics 
generated by the laser model. (B) Partial bifurcation diagram though the parameter continuation showing the curves for 
pitch- fork (PF) and Andronov-Hopf (AHq) bifurcations for the equilibrium state, O, and another similar supercritical one for 
Oi 5 2- The homoclinic curve, Horn begins from the codimension-two point, BT for the Khorozov-Takens bifurcation and ends 
up at the resonant saddle point. (B) Elevating a = 2 makes the Horn turned by a saddle point in the parameter plane and 
terminate at the primary T-point. 



the curve is continued, the homoclinic butterfly undergoes another codimension-two orbit-flip bifurcation 
so that the separatrices loops of the saddle, O become non-orientated. As the result, further down the curve, 
each loop gains an extra turn around the incoming separatrix of the opposite saddle-focus, i.e. becomes a 
double-pulsed one with the {1,-1,0} kneading. Depending on the parameter cut a there are two scenarios 
for termination of the curve, Horn, in the (a, b) diagram (compare the bifurcation diagrams in panels B 
and C of Fig. 17): first, when a — 1.5 it terminates at the codimension-two T-point corresponding to 
the heteroclinic connection between all saddle equilibria, O and as shown in panel B. The second 
scenario for is less predictable at a — 2: the branch, Horn, of double-pulsed separatrix loops ends up at 
the codimension-two point of the resonant saddle with the zero saddle- value (panel C). The answer to the 
question what makes the curve change its destination is a saddle point in the parameter diagram that the 
kneading scan reveals in Fig. [13] By varying the a-parameter cut in the 3D parameter space, this bifurcation 



curve is destined by the saddle to finish at either terminal point, see details in Shilnikov et al, 1993 . In 



the case where it spirals onto the T-point, there is another bifurcation curve corresponding to the same 
{1,-1,0} kneading that connects the codimension-two orbit-switch point and the point corresponding to 
the resonant saddle located on the curve Horn. 

The panel A in Fig. [IT] represents the kneading scans of the dynamics of the laser model which is 
mapped onto the (a, /3)-parameter plane with g = 50 and a — 1.5. The scan is done using the same 50 
kneading entries. It has the regions of chaotic dynamics clearly demarcated from the solid color windows 
of persistent kneadings corresponding to trivial attractors such as stable equilibria and periodic orbits. 
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The region of chaos has a vivid fractal spiral structure featuring a chain of T-points. Observe also a thin 
chaotic layer bounded away from the curve Horn by a curve of double-pulsed homoclinics with the kneading 
{1,-1,0} connecting the codimension-two points: the resonant saddle and the orbit-flip both on Horn. One 



feature of these points is the occurrence of the Lorenz attractor with one or more lacunae [ Afraimovich 
\et al\ [1983} Shilnikov, 1986, 1993; Shilnikov et a/., 1993] . Such a strange attractor with a single lacuna 
containing a flgure-8 periodic orbit in the phase space of the given laser model is shown in panel A of 
Fig[l6} 



7. Conclusions 

We have demonstrated the new proposed computational toolkit for thorough explorations of chaotic dy- 
namics in three exemplary models with the Lorenz attractor. The algorithmically easy though powerful 
toolkit in based on the scanning technique that maps the dynamics of the system onto the bi-parametric 
plane. The core of the approach is the evaluation of the kneading invariants for regularly or chaotically 
varying flip-flop patterns of a single trajectory - the separatrix of the saddle singularity in the system. In 
the theory, the approach allows two systems with the structurally unstable Lorenz attractors to be conju- 
gated with the mean of a single number - the kneading invariant. By using ready- for- use tools in Matlab, 
we could have the parameter plane of the model in question be foliated by the level curves of distinct colors 
corresponding to distinct values of the kneading invariants. The kneading scans revel unambiguously the 
key features in the Lorenz-like systems such as a plethora of underlying spiral structures around T-points, 
separating saddles in intrinsically fractal regions corresponding to complex chaotic dynamics. We point out 
that no other techniques, including approaches based on the Lyapunov exponents, can reveal the discovered 
parametric chaos with such stunning clarity and beauty. Figure [18] for the Shimizu-Morioka model shows 
that a fine scan based on the finite-time Lyapunov exponents is able to indicate some presence of spiral and 
saddle structures and differentiate between the chaotic dynamics due to the Lorenz attractor and those 
due to additional degrees of instability brought in by saddle-foci. 




1st Lyap.exp. 
1st Lyap.exp. >0 



1st Lyap.exp. <0 
2nd Lyap.exp. 



parameter a 

Figure 18. Fine scan based on the Lyapunov exponents indicating the presence of the saddle (in the Lorenz attractor 
region shown in cold blue) and spiral structures (in the reddish regions with larger Lyapunov exponents for wild chaos due to 
saddle-foci) in the (a, A)-parameter space of the Shimizu-Morioka model. 



The kneading based methods should be beneficial for detailed studies of other systems admitting rea- 
sonable symbolic descriptions. It bears an educational aspect as well: the kneading based scanning can 
be used for in-class presentation to reveal the fascinating complexity of low-order models in the cross- 
disciplinary field of nonlinear dynamics. The bi-parametric mapping technique can be easily adopted for 
parallel multicore GPU platforms allowing for ultra-fast simulations of models in questions. Additional 
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implementation of high-precision computations of long transients shall thoroughly reveal multi-layer com- 
plexity of self-similar fractal patterns near T-point vortices. In future research we intend to enhance and 



refine the toolkits for exploration of other symmetric and asymmetric IShilnikov & Shilnikov 



terns of differential and difference equations, like 3D Poincare mappings |Shilnikov et al. , 1993 



|Shilnikov 8z Shilnikov 


, 1991] sys- 


IShilnikov et al, 1993; 


Gonchenko| 



et al. , 2005 , including 4D models with saddle-foci, that require two and more kneading invariants for the 



comprehensive symbolic description. 
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